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TECHNICAL NOTE 4275 

DYNAMIC STABILITY OF VEHICLES TRAVERSING ASCENDING OR 
DESCENDING PATHS THROUGH THE ATMOSPHERE 
By Murray Tobak and H. Julian Allen 


SUMMARY 


An analysis is given of the oscillatory motions of vehicles which 
traverse ascending and descending paths through the atmosphere at high 
speed. The specific case of a skip path is examined in detail, and this 
leads to a form of solution for the oscillatory motion which should recur 
over any trajectory. The distinguishing feature of this form is the 
appearance of the Bessel rather than the trigonometric function as the 
characteristic mode of oscillation. 


INTRODUCTION 


The speed of many rocketcraft and some aircraft is 60 great that 
they possess ample kinetic energy to permit ejection from the atmosphere. 
With the satellite vehicle, the ballistic rocket, and the skip rocket, 
such trajectories are, of course, intentional. With high-speed aircraft 
and with glide rockets, although sojourns into space may purposely be 
made, they may also accidentally occur if, when in climbing flight, the 
dynamic pressure falls so rapidly that the vertical momentum cannot be 
checked to prevent such a course. 

When a vehicle is in space it may have arbitrary orientation with 
respect to the flight path and this orientation may change with time. In 
the absence of air, the history of this orientation can only be changed 
by the use of some form of reaction control. However, as the vehicle 
descends into the atmosphere it will develop motions in response to the 
aerodynamic forces and moments which come into play. If the vehicle 
possesses positive static and dynamic stability, these motions will sub- 
side into the types with which we have become accustomed from more con- 
ventional aircraft. The whole history of the motion is of interest to 
the designer, both for guidance system and structural design and in 
determining the feasibility of human habitation. 

The dynamic motion of ballistic vehicles has been the subject of 
several theoretical investigations. Oswald (ref. 1) considered the 
motions of an unguided rocket in its flight out of the atmosphere. The 
general problem of a ballistic missile entering the atmosphere was treated 
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by Friedrich and Dore in reference 2. The analysis of reference 3 treated 
the entry problem -with certain simplifying assumptions regarding the tra- 
jectory and the atmosphere in order to clarify the nature of the problem. 
For add the above analyses, the zero-lift trajectory was used. 

It is the purpose of this paper to investigate the motion history 
for the cases for which aerodynamic lift is important in changing the 
trajectory. In common with the analyses cited above, it will be assumed 
that the trajectory lies in a vertical plane and that disturbing forces 
are symmetric with respect to the vehicle's vertical plane of symmetry. 
This reduces the problem to one of longitudinal motions only. In order 
to illustrate the mechanics of the problem more clearly, the analysis is 
concerned first with the specific example of a vehicle which descends 
through the atmosphere on a skip trajectory, executes its turn, and is 
ejected again from the atmosphere. A detailed study of this problem 
reveals the nature of the essential parameters, and this leads to a more 
general form of solution which can be used to calculate the oscillatory 
history over any given trajectory. 


SYMBOIS 


A 

CD 

cl 

Cm 

C L, 


a 


"“a 


Cm d 


■'mq 


reference area 


drag coefficient. 


drag 

qA 

lift 

qA 


lift coefficient, 
pitching-moment coefficient, 

rate of change of lift coefficient with angle of attack, 
rate of change of pitching-moment coefficient with angle of 

/ay 


attack. 


\cScx 


s ) 

/a-*o 


rate of change of moment coefficient with time rate of change 
of angle- of -attack parameter 


ad ( dC m \ 


rate of change of moment coefficient with pitching velocity 


parameter % 



tilt-* 
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C™ rate of change of moment coefficient with control deflection. 

N -*>o 

e Naperian base 

g acceleration due to gravity 

I pitching moment of inertia about center of gravity 



Bessel functions of first kind of zero and first order 


K x 

Ke 


l 


m 

1 

R 


s 

t 

V 

¥ 


Yo( )l 
*i( >J 


a 


3 


dynamic stability parameter (eq.. (15)) 
static stability parameter (eq.. (2l) ) 
lift-drag ratio 

body length and reference length for moment coefficient 
evaluation 

vehicle mass 

dynamic pressure , ^ pV 2 

radius of earth 

distance measured along path of static trajectory 
time 

flight speed 

flight speed on entrance to the atmosphere 
vehicle weight 

distance measured along earth 
altitude 

Bessel functions of second kind of zero and first order 
angle of attack (sketch (a)) 

initial value of oscillatory angle of attack on entrance to 
the atmosphere 

density parameter (eq.. (8)) 
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7 flight path angle, a - 9 (sketch (a)) 

7 e static flight path angle on entrance to the atmosphere 

7 S static flight path angle 

8 control deflection angle 

9 angle of pitch (sketch (a)) 

| dummy variable 

p air density 

p Q air density at sea level 

a radius of gyration . . .. 

w oscillation circular frequency, radians/sec 


ANALYSIS 

The plan of the analysis to be undertaken is as follows: We shall 

present first the equations which apply generally to the longitudinal 
motions of a vehicle moving through the atmosphere at high speed. Then, 
several sections will be devoted to a study of the equations as they 
apply specifically to the skip trajectory, as this will reveal a number 
of important features and permissible simplifications. In a final sec- 
tion, these features and simplifications will be exploited to derive a 
general form of solution for the oscillatory motion which is applicable 
over any given trajectory. 


Equations of Motion 


\ 

X 



Consider a vehicle moving 
through the atmosphere at high 
speed, which at time t is 
approaching the earth along a 
flight path depressed from the 
horizontal by the positive angle 
7 (see sketch (a)). The full set 
of equations of motion defining 
the vehicle 1 s path and the oscil- 
lations of the vehicle about that 
path may be written as 


Sketch (a) 
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-mV - CpqA + mg sin 7 = 0 
mVy +0 L qk+m(^~ - g)cos 7 = 0 


I' 


IS - qkl > C m = 0 


CD 


wherein 


6 = a- 7 

Cl = C^a 

^Cm = C mcL a+ ~y~ + Cj^ -y- + CjQgS 


and a dot over a symbol denotes differentiation with respect to time. 
Because of the dependence on time of density p and velocity V, equa- 
tions (l) are nonlinear; their solution cannot be effected analytically 
unless simplifying assumptions are made. Herein we assume the following: 

1. The aerodynamic coefficients in equations (l) are essentially 
independent of Mach number, a Justifiable assumption for the high flight 
velocities of interest here. 

2. The drag coefficient Cp is independent of angle of attack and 
pitching velocity, a justifiable assumption so long as these variables 
remain small. 


3. The equations may be broken down into two sets of equations, 
one set defining the “static trajectory” of the vehicle's center of 
gravity, the other describing the osc i n ations of the vehicle about its 
center of gravity as well as the oscillatory motions of the center of 
gravity about the static trajectory. The justification for this proce- 
dure will be evident from the following considerations. Each of the 
dependent variables a, 7, and 9 in equations (l) is comprised of a 
uniform variation with time plus an oscillatory variation. Hence, let 
us make the following substitutions: 


7(t) = 7 S (D + 7 0 (t) 

a-(t) = Os(t) + cc^Ct) > 

e(t) = 0 s (t) + e D (t) 


( 2 ) 


wherein the subscript s denotes the uniform part of each motion and o 
the oscillatory part. Now, under the restriction that y Q be a small 
quantity, so that 
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« 1 


sin y g [ 1 + 


7 - 0 ), 

yj 


sin / 


cos 7 1 + — « cos 7 


equations (l) may be grouped as follows: 


-mV - CjjqA + mg sin 7 = 0 


[mV7 s + qACp^ag + m(^$y - g)cos 7J + |mV7 0 + qAC La a 0 J = 0 
[l0's “ lAZ (Cjj^cXg + Cjjj^ -y— + Cjjj^ -y- + CmgSylJ + 

[l0 O -^(c^Oo + C^ ^)j = 0 


> ( 3 ) 


In each of the last two equations, the second bracket is a sum of oscil- 
latory quantities whereas the first bracket is not; under the assumption 
that over a given cycle of oscillation the velocity V and dynamic pres- 
sure q remain essentially constant, the average of the second group of 
terms over the cycle must be near zero. Hence, so must the average of 
the first. But since the behavior of the first group of terms is not 
oscillatory, in order that its average be near zero its value at any 
time must be near zero. Therefore, we can assume that each group of 
terms in equations (3) is independently zero. We have then for the non- 
oscillatory or ” static” trajectory _ 

-mV-C^qA+mg sin 7 g = 0 1 


mV7 B +qAC La o s +m(^ ; - g)cos 7 S = 0 
I0 8 -qAZ^as + C^ & -f- + ^ = 0 


> (4) 


Qq - ^3 ~ 7 B J 


and for the oscillatory motion 
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mV7 0 + qACL^cto = 0 

.. f 0 O Z C ir>T>\ 

lfe 0 ” <3^1 ^Cjjj^iXq + ~y i- ~~Y~) ~ ® ^ 

9 o = a o " 7 0 


( 5 ) 


The set of equations (5) can he combined into a single equation for the 
osci ll atory angle of attack of the form 


wherein 


cto + f i(t )cio + — ® 


f i( t ) = °La _ ( C m cl +C mi) ^§1 

_d A pVA\ _ Cm q C 4q /pffiA 2 

dt y/La 2 m/ Im 2 J 


( 6 ) 


f 2 (t ) 


pV^AI 
^ 21 


> (7) 


Skip Trajectory 


The differential equations (4) to (6) are applicable to the analysis 
of arbitrary static flight paths and their accompanying oscillatory his- 
tories. We now consider the application of equations (4) to (6) to the 
specific case of a skip trajectory under the simplifying condition that 
L/D is constant. 

Static trajectory .- Since Cp has been assumed constant, the 
specification that L/D be a constant means that % is also fixed. 

This permits the variation of flight path angle with time to be deter- 
mined from the first two of equations (4). The third equation specifies 
how the control motion is to be programmed in order to maintain constant 
L/D. We shall not be concerned with the latter problem here. 


The static skip trajectory for constant lift-drag ratio ha6 been 
determined in reference 4 under the assumption that: (l) the air density 

varies with altitude in the manner 


P = 



(8) 


wherein p Q and p are constants, and (2) the effect on the path of 
gravity and the term V 2 /(R+y) may be ignored. The first assumption is 
an adequate one (see ref. 5) over the range of altitudes for which the 
aerodynamic forces are important in determining the path, provided the 
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constants p Q and p are appropriately chosen. In most cases the values 

* 

p Q = 0.003^ slugs /ft 3 
p 22,000 ft 

are satisfactory and -will he used in this report. The second assumption 
is necessary to keep the functional relations single. It Is a satisfac- 
tory one when the flight velocity is high and the aerodynamic forces 
experienced during the skip are large (see ref. 4 ). For present purposes 
It Is felt to he sufficient, since our concern is more with the nature of 
the motion problem rather than with an accurate evaluation of the motion. 
Nevertheless, the adequacy of this assumption will he reviewed later in 
this report as it affects the quantitative results. 


If one defines y s as the instantaneous static flight path angle 
and y E as the static flight path angle at the "entrance" to the atmos- 
phere, then the analysis of reference 4 gives the velocity and density 
dependence on y s as 


y s -?e 
l/d 

V = V E e 

-py 2pm / v 

P = Po e = ( GGS 7 S “ GGS 7 e) 

Cl = CLqOs = constant 


( 10 ) 


The coordinates of the flight path can he 



In - 


c LPo A 


2pm(cos 7 s - cos 7e ) 




found from equations 

-\ 


(10 ) to he 


x = 



I 1 / tan y -%/2 + tan[ y s /2 [\ 

7s ' + tan 7 E Zn \^ban 7e /2 - tanj 7s /2 \J 




J 


(H) 


In the equation for x the minus sign is used for the incoming path and 
the plus sign for the outgoing. 


The path length 
shown to he 


s measured from the bottom of the skip can he 


s = =F 


In: 


(tan 7 e / 2 + tan 


P sin 7 E \tan 7 e /2 - tan 


7 s /2 


7 b /2 


( 12 ) 


wherein, as with x, the minus sign is used for the incoming path and 
the plus sign for the outgoing. 
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Oscillatory motion .- Having eliminated a s from equations (10) a nd 
7 0 from equations ( 5 ), henceforth we shall omit the subscripts s and o 
as no confusion should ensue. It will be understood that in referri ng 
to a we mean the oscillatory angle of attack, whereas in referring to 
7 we mean the flight path angle of the static trajectory. 


Consider now equation ( 6 ) for the oscillatory angle of attack, 
will be found convenient to consider the independent variable as 7 
rather than t. Using the relations 


da, 

dt 


a, 


(7) & 
dt 




It 

(= 7 b ) 



d7 

dt 


7-7e 

l/d 

"PV E e 


(cos 7 - cos 


7e) 




(13) 



3ln 7 ^ 

^dt/ \L/D cos 7 - cos 7 -g/ 


J 

we find, after some manipulation, that equation ( 6 ) becomes 


a"(7) +f 3 (7)0' , (7) +f 4 C7)ct'(7) = 0 


(14) 


where 


f 3(7 ) 
f4(7) 


Ki 


sin 7 
pxr _ _ _ / 

1 cos 7 - cos 7 g 

1 [52 , ^ _1_ AY] 

(cos 7 - cos 7 g) L^L s n 7 ^ (32 \o J J 

r C lq, 1 QmqpLq. AY' 

L°l L/D Cl 2 W . 

1 r_l_ C ^g + C ^a AY' 

2 Ll/d " C L + Cl [aj _ 


> (15) 


J 


and the primes indicate differentiation with respect to 7 . Fi nally , 
following the method of Friedrich and Dore (ref. 2) , we remove the first- 
derivative term in equation (l4) by means of the transformation 


q(7) = a(7)e 



(16) 
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whereupon equation (l4) becomes 

tf’(r) +M(y)B(y) = 0 ( 17 ) 

with 

M(r) = f 4 (r) - \ f 3 '(r) - | (18) 

Alternate solutions of equation (l4) are presented in the following two 
sections. 

Method of Friedrich and Dore 


Before proceeding to the solution of equation (l4) which was used in 
the body of the report, it is useful to present an alternate solution 
based on the method of Friedrich and Dore (ref. 2). Such a solution is 
applicable under somewhat more general circumstances than those of prin- 
cipal concern here and, hence, may prove valuable in other applications. 


In the method of reference 2, it is necessary to assume that some 
technique is available for calculating the motion a,(y ) from the entrance 
condition down to some point y p at which point a(y p ) is at the peak 
of a cycle. Let this value of a be called' cop. Then an excellent 
approximate solution for the envelope curve of the a oscillation from 
this point forward is 


a max(7) - Op j 


[M(7 P )1 

1 

4/bos 7 p - cos 

lM(7) J 

^ cos 7 - cos 


v \ 2 

4 ) 


(19) 


where M(y) is given by equation (l 8 ). Equation (19) can be used to 
est im ate the magnitude of the a motion under fairly wide variations of 
the atmospheric, aerodynamic, and inertial parameters. Further, using 
equation ( 19 ) and following the arguments of. reference 2 , one can go on 
to compute the detailed history of the angular oscillations , frequencies, 
etc. These results will not be presented here in full, since their deri- 
vation is well described in reference 2. Also, primarily because the 
expression for M(y) is very complicated, the full expressions are some- 
what unwieldy. For the conditions of concern in this paper, however, it 
will be shown later that M(y) is given with good accuracy for all values 
of y, except those very near y E , by 


where 


M(r) 


Kg 

cos y - cos /g 



( 20 ) 


( 21 ) 



MA.CA. TBT 4275 


11 


Since in equation ( 19 ) it must be assumed that 7 $ 7 < 7 , the approxi- 
mation (20) can be used, and equation (19) simplifies uo 


° Wi.y ( 7 ) — 



( 22 ) 


Alternatively, the form of the dependency of a ^.y on atmospheric density 
can now be demonstrated explicitly by use of equation (10); thus 


rp(7 P )i4 

°ffla x(7) = op e 


Kx(7p-r) 


( 23 ) 


Finally, under the assumption that the static flight path angle does not 
change radically from its initial value, applying the small-angle approxi- 
mation to equation (10) permits equation ( 23 ) to be ■written as 


omax(7 ) = 


1 


^[H7rr e,Vl[p ' 7, ' p<rp)1 ’ r>0 

^J”p(7p)_J* e ^iC P ( 0 ) “ p(7p)] e Ai[p(0) - p(y) ] 


7 <0 (24) 


• with 


*1 = 


AC^Kj. 
2pm sin 7 E 


(25) 


The first of equations (24) applies to the descending phase of the tra- 
jectory, the second, to the ascending phase. To the same order of accuracy 
as equation (22), the variation of oscillation frequency (in radi ans /sec) 
is found to be 


w = PV E jK^cos 7 -cos 7 e ) 


7~7 e 

l/d 


( 26 ) 


Alternatively, using equations (10) again, we may put equation ( 26 ) in 
the form 


“ ^( 7 ) | v/pTtT ( 27 ) 

The variation of maximum angular velocity may be found approximately from 
equations (24) and ( 27 ) according to 

ct max(7) = a max(7) w (7 ) 


( 28 ) 
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It will be noted in equations ( 24 ), equation ( 27 ), and so, also in 
equation (28) , that in effect the details of the particular trajectory 
under study have been eliminated. Hence, these expressions may be expected 
to yield at least qualitative estimates of cy^y , dn fly , and oj for other 
trajectories not too dissimilar from the one considered here. This point 
will be exploited later as a means of studying the effects of gravity on 
the motion. 


Bessel Function Solutions 


We now develop an approximate solution to equation ( 17 ) from whose 
form the main characteristics of the motion can be deduced with relative 
ease. Basically, the solution depends on finding an adequate approxima- 
tion to M(7) over the entire range of 7. This can be done by consider- 
ing the relative magnitudes of the terms in M(y ) for some representative 
values of the parameters. Such a study was carried out (see appendix) 
and it was found that for all values of 7, M(7) can be approximated with 
excellent accuracy by 


M(7) 


K 2 I - cos 7 cos 7 E 

cos 7 - cos 7g 4(cos 7 - cos 7g) 2 


( 29 ) 


Wow since KA is of the order of 10 s , it is clear that the second term 
in equation (29) is trivial compared to the first for all values of 7 
except those very near 7^. Since the Friedrich-Dore method may be 
applied only after the motion has completed at least one half-cycle, 
values of 7 near y E are automatically excluded; hence the second term 
may be discarded and we have the approximation used previously in equa- 
tions ( 20 ) and ( 22 ). Here, however, we seek another approximation which 
will include the first half-cycle. Such an approximation can be made if 
we require that 7^, the entrance angle, be not too large. In this case, 
using the small -angle approximation, we have for 7g>7>0 


cos 7 - cos yg = 2 sin 




( 30 ) 


Likewise, in the range 0 > 7 > -y E , the appropriate approximation is 

cos 7 - cos 7 E « (7 E +7)sin 7 E ( 31 ) 

As for the second term in equation ( 29 ), since it is of significant 
magnitude only for values of 7 near 7^, a valid approximation to the 
numerator is required only for those values of 7. Hence, 

1 - cos 7 cos a 1 - cos 2 7g = sin 2 7g 


( 32 ) 
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Finally, therefore, using equations (30), (31) , and (32) , we have as 
approximations to M( 7 ): 


7 E >7 >0 i 



> 


M(y) 


+ 1 

0 >7>-7 e ; 

7 e “7 

M7 e “7) 2 


M ( 7 ) 

*2 

, 1 


7 E +7 

+7) 2 

where 



k 2 


k 2 “ sin 7 e 


(33) 


The first of equations (33) is to he used for that part of the trajectory 
in which the vehicle descends through the atmosphere j the second, after 
the vehicle "begins to ascend. The two phases will now he considered in 
turn. 


Descent phase .- Letting the variable be rj = 7 E ~ 7 , and substituting 
the first of equations (33) in (17 ) } we have for the equation of motion 
in the range 0 <tj <7g, 

+(^ + i)a(n) = 0 (3*) 

The solution to this equation is known (cf . , ref. 6) and gives 



where <J 0 (z) and Y 0 (z) are, in Watson’s notation, the zero-order Bessel 
functions of first and second kind, respectively. The exponent in equa- 
tion (l6) may be evaluated by straightforward integration, whereupon, 
reverting to the original variable, we have for the a motion 

a(y) = e Kl ^ 7 E 7 ^a 1 J Q (2 >v /K 2 (7 E - 7 /) + ( 7e -/))J (36) 

The constants a x and a 2 are of course to be evaluated by specifying 
the initial conditions. If we specify that the vehicle enter the 
atmosphere at initial angle of attack Og and without angular velocity 
then the second term in equation (36) vanishes; hence the a motion in 
the descent phase is given by 
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*(y) = c£Ee Kl ^ 7E 7 ^Jo(2 Jk 2 (7 e -7)) i 


7 e >7 >0 


(37) 


Ascent phase .- With the variable now T] = / E + 7, substitution of the 
second, of equations (33) in (17) shows that the equation of motion is of 
the same form as equation (34). Hence, the solution for S(ti) is again 
given by equation (35) • After evaluating the exponent in equation (l 6) 
(with the limits of integration, 0> | > 7) we get 


a(y) = e Kl7 cJ 0 ^jK 2 (y E + y)^ + dY 0 (2 N /K 2 (7 E +7)) 


(38) 


Ihe constants of integration again may be evaluated from the initial 
conditions, which are, of course, the end values of the descent phase. 
From equation (37); with y = 0, we have 


a(0) = a E e K:L7 Ej 0 (|i) 

a'(0) = a E e Kl7 E^-K 1 J 0 (p) + Ji(n)] > 

P = 2/Ki y^ 


(39) 


Using equations (39) together with (38), we may evaluate the constants 
c and d, whereupon the a motion in the ascent phase is found to be 


a (y) = 0^ 


Ki(7 B -7) 


»rrJo'(n)Ti(n) + Ji(n)Yo'Cn)l T L r—r— x\ 

lfe(7)V*y ■ - 5aiK 2W, < + r V 


f; 




Jo(ji)Tx(ft)-Ji(ui)Yo(n)J 


3 (2 ^s(7 B +7))} i 


0 >7>“7 e 


m 


Finally, it is clear that the same technique used to cast equa- 
tions (24) to (28) solely in terms of atmospheric density may be U6ed 
in the results (37) and (40) as well. Thus , we have as an alternate 
solution for the motion in the descent phase 


0,(7) = Oge 


A-,p(7) 


'0 ^n/^ 2P (7 )^ > 


7%>7>0 


(41) 


Likewise, the solution in the ascent phase becomes 
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_ >i[p(0)-p(7)]> 1 p(0)rrJ o (n)Y : 

a(r) ■ ^ e te* 




where 



0>7>-7 E 


(42) 


(43) 


DISCUSSION 
Special Cases 


It is instructive to study equations (37) and (40) under some 
special conditions , as the character of the motion then can be deduced 
almost by inspection. Thus, let us assume that the dangling coefficient 
Ki in equations (37) and (4o) is small enough to be neglected, and let 
us examine in turn the three special cases; (l) J x (p) = 0; (2) J 0 (i-l) = 0 
(3) Jo(p)Yi(p)+Ji(p)Yo(p) = 0. 

Case (l); Jx(p) = 0 « ~ Since in the descent phase a( 7 ) is dependent 
only on the Bessel function J 0 (2 " y) ) > and since Jx(n) = -Jo'fp) 

the condition Ji(p.) =0 is to be interpreted as meaning that a is at 
the peak of a cycle at the instant the static trajectory is at its lowest 
point. The a motion in the ascent phase then reduces to 


a(7) = ogJo^Ka (7 E + 7 )^) ; 0>7>-7 E 


(44) 


Hence, we find that the motion in the ascent phase is simply the mirror 
image of the motion in the descent phase; that is a ( 7 ) is an even 
function about 7 = 0 as shown in sketch (b). 

Case (2); J Q (p. ) = 0.- With J 0 (p) = 0, a passes through zero at the 
lowest point of the static trajectory. Equation (40) reduces to 
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o o(r) = -c^J q (2n/k2(7e + 7)) i 0>7>-7 E 


( 45 ) 


Hence, the motion in the ascent phase is the negative mirror image of the 
descent phase j that is, 0 .( 7 ) is an odd function about 7 = 0 as shown in 
sketch (c). 




Case (3); Jp(n)Yi((i) + Ji(m.)Y 0 (m.) = Q .- Examination of this Bessel 
function combination for large values of the argument reveals that it is 
zero when a has completed either one-quarter or three-quarters of a 
cycle at 7 = 0. Equation (4o) reduces to 


a 


= 




0>7> -7 E 




Note that, unlike the two previous cases, in this case the vehicle leaves 
the atmosphere with an angular velocity, as shown in sketches (d) and (e) , 
since the Y 0 Bessel function approaches infinity as the argument 
approaches zero. Although a , 1 ( 7 ) as - 7 E , nevertheless, the 



Sketch (d) 


Sketch (e) 
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angular velocity a(t) is finite, since y(t)-*0 as y-^-y-g. By the 

appropriate l imi ting process, the value of d(t) acquired by the vehicle 
as it leaves the atmosphere may be calculated to be (for Ki » 0 and 
case (3)) 

- 2 7 E 

(S)exit 6 ^ a E PY E sin 7 E <*) 


More generally, when Ki is not neglected, and for an arbitrary value 
of a(o), the angular velocity at exit is 


da'N 

dty exit 


-2y 

Bin 7E J oC^)Ji(^)e 



m 


Hence, for all but the special cases (l) and (2) we may say that the 
vehicle will tumble as it leaves the atmosphere, with a rate given by 
equation (4-8) . 


Oscillation History for Representative Vehicle 


Let us now use the full equations (37) and (4-0) to determine the 
oscillation history of a vehicle representative of a type which might be 
suitable for flight at very high velocities and altitudes. The airframe 
selected consists of a low aspect ratio triangular wing mounted on a 
conical body, 50 feet in length. We assume it capable of developing an 
L/D of 6 at angle of attach of 0.1 radian. Calculations indicate that 
a representative value of the damping factor Ki for such an airframe 
is of the order of -10; the static stability factor Kg is of the order 
of 3x10 3 for a static margin of 3 percent. Let us also specify that the 
vehicle have an over-all range of 4-000 nautical mil es, so that, using 
results derivable from reference 4-, we have that the vehicle enters the 
atmosphere on the skip phase of its trajectory with an initial velocity 
of 14,500 feet per second and with the flight path inclined initially at 
y E = 12.2°. 


In figure 1 is shown the history of the angle-of -attack oscillation 
together with the static skip trajectory. Also shown dotted is the his- 
tory with the damping factor Ki taken an zero. It will be noted that 
on the descent phase, even with zero damping, the motion subsides very 
quickly because of the powerful influence of increasing atmospheric 
density on the aerodynamic restoring moment. On the ascent phase, how- 
ever, the aerodynamic damping is surprisingly effective In keeping the 
amplitudes small as the vehicle gains in altitude. Also, it should be 
noted that the value of |_l for this example (eq. (43)) is such that cc 
passes neither through zero nor a peak at the bottom of the skip; hence , 
the vehicle will tumble as it leaves the atmosphere. The tumbling 
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rate does not appear to "be very serious, however; for the present example, 
with og = 5° and even with Ki taken as zero, it is calculated to be of 
the order of one revolution every thirty minutes. With the damping factor » 

included, this result is diminished seventyfold - to the order of one 
revolution every thirty-five hours. 


General Form of Solution 


If one examines the results of references 2, 3> and those of this 
paper (all dealing with the oscillatory motions developed by vehicles on 
entering the atmosphere), one finds a feature in common that is a depar- 
ture from the familiar characteristics of conventional aircraft motions; 
namely, that the characteristic mode of oscillation is described by a 
Bessel function rather than a trigonometric function. Physically, the 
reason for this is clear since, as a vehicle descends through the atmos- 
phere, the rapidly increasing atmospheric density causes the aerodynamic 
restoring moment to act in the manner of a stiffening spring. Thus, the 
oscillation must tighten both in amplitude and frequency, and this is 
precisely the behavior described by the Bessel function. 

It already has been demonstrated (eqs. (4l) and (42)) that the Bessel 
function solution for the oscillatory motion during the skip trajectory 
could be put in a form that is, in effect, independent of the type of 
trajectory. We wish now to show that this form will in fact recur for 
any trajectory, or any part of a trajectory, in which the direction of 
flight does not change rapidly with altitude. To this end, let us con- 
sider equation (6) again, and now change the independent variable to 
altitude y. The result Is 

a"(y) +a‘(y) 7 »(y)cot 7 + ^ § “ y" s ' ln y f * y 2 S i n a y f s(y) a (y) 
with 

* 

fs(y) - - ir- 

and where we have retained only the dominant term In f 2 (y). Wow, in 
many cases (for example, on entering or leaving the atmosphere) y will 
change only very slowly with altitude. Hence, let us assume 

7 ’(y) * o 


= o 

m 


(50) 


) 


sin y * constant 


(51) 
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Equation (49) becomes 

a"(y) + a' (y)^ ^ + 20A 3 p(y) J + p 2 ?\ 2 p(y)a(y) = 0 

with 


p(y) = P n e 


-py 


A. = 


' 3 4pm sin 7 

-A 


[ _CL a + (a) (% +Cm d)_ 

AC L Ks 


2 p 2 mZ sin 2 7 CO Cma Oo ” “ " v ' 2 


2pm sin 2 y 


(52) 


(53) 


Note that X 3 and ~h 2 have the same form as ~K X and 1 2 , respectively, in 
equation (43) except that new A 3 does not contain the drag coefficient 
pi.nd y is not necessarily the entrance flight path angle. Making the 
same transformation as before , 


l 

"2 

a(y) = a(y)e 


J 7 x l 


^fjT + 2pA 3 p(|) 


4g 


(54) 


we have 


a” (y) +N(y)a(y) = 0 


(55) 


with 


N(y) = P 2 A 2 p(y) 


2pA 3 p(y)] - i 


i g + 2pA 3 p(y) 


(56) 


Now since A 2 is proportional to K 2 whereas the bracketed term is 
essent ially proportional to Ki, it is easy to show (see appendix) that 
the contributions of the two terms containing \ 3 and V are small 
compared with that of the first term. Retaining only the first tern in 
equation ( 56 ), we get for the differential equation 

a"(y) +P 2 7\ 2 p(y)s(y) = 0 (57) 

in which form the nature of the aerodynamic spring is clearly evident. 
The solution to equation (57) is 


a,(y) = c 1 J 0 (2^ 2 p(y))+c 2 Y 0 (2^ 2 p(y)) 


(58) 
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Consider now the exponential factor in equation (5^). Recognizing that 
— jjj-jj- can he written as ^ (in V), ¥e see that this term is immediately 

integrable. The exponential term becomes 



(59) 


Note that the damping factor 7\ 3 has been kept inside the integral, 
since by this means we extend the validity of the analysis to cases 
wherein 7\ 3 may be dependent on altitude. Then, combining equations (58) 
and (59 ) > we have as a general form for the oscillatory motion 

r y 

-0 J A 3 p(S)d| 
e 



CiJ 0 (WA 2 p(y))+ C 2 Y 0 (k/A 2 p(y)) (60) 


Equation (60) has the following virtues which make it particularly 
suited to studying the oscillatory history. First, it is clear that the 
effects on the motion of gravity, variations in the drag coefficient, 
and use of thrust are all automatically incorporated by using in equa- 
tion (60) values of V(y) derived from an accurate solution of the static 
trajectory equations. Second, one can account for the effects of varia- 
tions with altitude of the damping parameter ?\ 3 . Finally, equation (60) 
can be used to derive a simple criterion whose satisfaction ensures that 
the motion is convergent at any altitude. Thus, noting that the envelope 
of the oscillations varies according to 

-0 ap(£)d| 

a max(y) = G ^ S T“ ( 6l ) 

[q(y)l* 

we have that the motion is convergent during a descent provided that 
a , max (y)>0. This is ensured if 

-0A 3 (y)p(y) - ^ ' ^y^ ' >0 ( 6s ) 

Again, it should be noted that the relation (62) correctly includes 
gravity effects when these are included in q(y). In this respect, 
equation (62) is a more precise statement of a similar result presented 
by Allen in reference 3, in which gravity effects are neglected. 

Next, it is in order to show that the form of equation (60) does in 
fact reduce to the form derived previously (eqs. (4l) and (^2)) when 
gravity effects are not excessive. To see this, consider the first of 
the static trajectory equations (eqs. (4)). We have 
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i av QpA. 

V dy 2m sin y 



_ Ap(y) r _ W Sin y~ 
2 m sin 7 A q. 


( 63 ) 


From the second of these forms, it is clear that the gravity term can he 
ignored over any part of a trajectory in which the inequality 


_J1 » JL 

[sin 7 QdA 


m 


holds. Alternatively, as a more general expedient one may average the 
"bracketed quantity in equation (63) in some sense and call it an effec- ' 
tive constant drag coefficient, say Cpi e . In either of these circum- 
stances, equation (63} can he integrated. Using the latter, we get 


V = V(y x )e 


ACp e 

23m sin 7 


Cp(y) - p(yx)] 


(65) 


Combining this result with equation (60) for the case X 3 = constant 
then gives 


o-(y) = e AlP ^ CiJoj&y/AaP (y)) + C 2 Y 0 (k/?V2p(y)) 


(66) 


with 

= 4 pm sin 7 " pE tt + © (>1 + Si)] 


Clearly, equations ( 4 l) and ( 42 ) are in accord with this form merely with 
the substitution of y for 7 as the independent variable. Further, by 
making use of the properties of Bessel functions for large values of the 
argument, one can show that the amplitude and frequency of the oscilla- 
tion vary essentially as 


cc max(y ) — 


_ 9hl 


Xip(y) 


[p(y) l 4 


«(y) = m v ( y ^p( y ) 


in confirmation of equations ( 24 ) and (27). Finally, it should be 
emphasized that should the inequality ( 64 ) become an equality, Cp e 
to be put equal to zero in equation (66). This of course is also 
reflected in equation (60) by the fact that V becomes constant. 


( 67 ) 


is 
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Arbitrary trajectories .- Although equations ( 60 ) and ( 66 ) are 
strictly applicable only while y remains essentially constant, in many 
instances the fact that y changes will be of no consequence, since, if 
the exponent is negative and the density increasing, the motion damps so 
rapidly that by the time y has changed significantly a is only a small 
fraction of its initial value. This is the case for the ballistic trajec- 
tory (cf., ref. 3) • On the other hand, should the sign of y change over 
the trajectory, as it does for the skip path, or if greater accuracy is 
desired, then the trajectory can be broken into an appropriate number of 
straight-line segments (none having zero slope, however), over each of 
which equations (6o) and (66) apply strictly. The parameters Ai, A 2 , 

A 3 and constants Ci, C 2 will then change from segment to segment, the 
first three being dependent on the value of sin 7 (see eqs. ( 53 ) and 
( 66 )), the second two being determined so that the initial conditions at 
the start of a segment match the end conditions of the preceding segment. 
It is interesting to note that the results previously obtained as equa- 
tions (4l) and (42) for the skip 
trajectory can be reproduced by 
this procedure If, after deriving 
the equations corresponding to a 
static trajectory composed of two 
line segments inclined at 7 E and 
- 7 g, respectively, one then assigns 
values to the ordinate y corre- 
sponding to those of the actual 
trajectory. Sketch (f) illustrates 
the procedure, which obviously can 
be extended to apply to arbitrary 
trajectories . 

Effect of gravity .- The obvious advantage - of equations ( 60 ) and ( 66 ) 
is that they may be used with any static trajectory whenever the altitude 
and velocity histories are given. As an immediate consequence, we may 
now evaluate the errors in our former results (eqs. (37) and (4-0)) caused 
by having neglected the gravity terms in the solution of equations (4) 
for the static skip trajectory. 

In order to carry out such a study, the static trajectory equa- 
tions (4) with gravity terms included were solved numerically for the 
same conditions as were used for figure 1 . 1 Having the new altitude and 
velocity histories, we could then use equations ( 6 l) and ( 67 ) to estimate 
the variations in amnv and to corresponding to this trajectory. Results 
for the static trajectory and velocity histories with gravity included 
and omitted are shown in figures 2 (a) and 2 (b), and the corresponding 

iThe term mV 2 cos 7 s /(R+y) in the second of equations (4) was neg- 
lected in these calculations in order to isolate the effects of gravity. 
For velocities of the magnitude considered here, its effect on the tra- 
jectory is small compared with that of gravity, and in the direction of 
reducing the differences between the results with gravity terms Included 
and omitted. 





Q 


NA.CA TN 4275 


23 


variations of O mpy and to in figures 2(c) and 2(d). It will be noted in 
figure 2(a) that because of the initial gain in flight speed (fig. 2(b)), 
the true static trajectory descends sooner into the atmosphere and also 
deeper, by about 15*000 feet. The effects of these changes on the oscil- 
latory angle-of-attack envelope (fig. 2(c)) are evidently not excessive; 
however, the frequencies are significantly increased (fig. 2(d)). This 
is a reflection of the fact that whereas the envelope curve varies as 
the fourth root of density, the frequency varies only as the square root. 

We should note that the above effects of gravity diminish as the 
trajectory becomes steeper. To illustrate this, a set of curves paral- 
leling figure 2 has been prepared for the case of a skip rocket having 
an L/D of 3 rather than 6. In order to achieve the same 4000 nautical 
mile range, the entrance angle in this event must be increased from 12.2° 
to 18.6° and the entrance velocity from 14,500 to 18,200 ft/sec. The 
results are shown in figures 3 (a-) to 3(d), and it is clear that here the 
influence of gravity is not significant either as it affects the static 
trajectory (figs. 3(a) and 3 (b)) or the oscillatory quantities (figB. 3(c) 
and 3(d)). On the other hand, however, since the accelerative loads 
become very large during steep skips, it would appear that a shallow 
trajectory of the type shown in figure 2(a) might be a more likely pros- 
pect for practical consideration than the one shown in figure 3(a). For 
such shallow flight paths, in order to obtain a static trajectory of 
sufficient accuracy, the gravity terms probably ought not be discarded 
in equations (4) . Once having obtained the static trajectory, however 
(by whatever means), equation (6o) or (66) can always be used to 
compute the oscillatory motion. 


CONCLUDING REMARKS 


An analysis has been carried out of the oscillatory motions of a 
vehicle which traverses a path through the atmosphere at high speed. 

The main results may be summarized briefly as follows: For any part of 

a static trajectory in which the flight path angle does not vary greatly 
with altitude, the oscillatory angle of attack and circular frequency 
will, have the form 


a 


(y) = 


-P f 7 ' Asp(i)dg 


•My) 


[c i Jo(2 W) + C 2 Y G (2 N /X 2 p(y)^ 


co(y) = KV(yWp(y) 


where p(y) and V(y) are, respectively, the atmospheric density and 
flight speed as functions of altitude. Since a trajectory can always be 
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decomposed into straight-line segments, over each of which the above 
equations apply, it is possible to say that given the altitude and veloc- 
ity history of a static trajectory, the corresponding oscillatory motion 
can always be computed by the use of these equations. 


Ames Aeronautical Laboratory 

National Advisory Committee 
Moffett Field, Calif., 


for Aeronautics 
Feb. 25, 1958 
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APPENDIX 

SIMPLIFIED EXPRESSION FOR M(r) 


The purpose in this appendix is to demonstrate that, for representa- 
tive values of the parameters, the coefficient M(/) can in fact he 
simplified to the form given in the text as equation ( 29 ). 

Written in full, Mfr) is expressible as 


M(7 > - 7 ? ) t 0111 7 (s 2 + Kl ) + I 008 7 +K2 ] + 

1 - cob y cos y E _ [% _i A\* 2 1 

4(cos 7 - cos 7 e ) 2 LQL L/D Cl 2 VV 1 J 


(Al) 


where 


K 


l 


1 r 1 + 

2 [l7d Cl 





_1_ (if 

c L pz W 


Hence, (comparing equation (Al) with equation (29)) we must demon- 
strate that the quantity within the brackets of the first term in equa- 
tion (Al) is essentially and that the third term is negligible 

compared with the other two. 


Consider now an airplane configuration which might be selected for 
flight at very high altitudes and velocities. Let it have a low aspect 
ratio triangular wing mounted on a conical body. A generous L/D for 
such a configuration is 6. As further specifications, let the vehicle 
have an over-all length of 50 feet, a trim angle of attack of 5*73° and 
(to prejudice the case against the argument) a static margin of only 
1 percent, measured in body lengths. Then 

*-« 1 


CL 


10 


_ _1_ 
cl " “ io 




(A2) 
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and, using linearized potential theory for triangular wings having 
supersonic edges, 

^ +C n& 1 

Cl ~ " 2 




cl 


(A3) 


Assuming the body to be uniformly solid, we find the moment of inertia 
about the center of gravity is essentially 


I as gg ml 2 


so that 


e-j 


« 22 


(a 4 ) 

(A5) 


Hence, with 


-L _ 22,000 = 44Q-- 

pZ " 50 440 




Kg = 968 

K x = -10.4 


CLq 1 

CL L/D 


+ 


CmqCLo, f 

1TW 


ffj + K 1 


> 


- 22.1 




j 


(A6) 


Obviously, then, even with an extremely small static margin, the influence 
of K 2 is overriding in the numerator of the first term of equation (Al) . 
Further, since the smallest the first term can become is Kq/( 2 . - cos 7g), 
the third term is negligible compared with the first for all values of 7. 
Therefore, M(7) is given with negligible loss in accuracy by 


M( 7 ) 


1 - cos 7 cos 7 e 

cos 7 - cos 7g 4 (cos 7 - cos 7g) 2 


(A7) 
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Range from lowest skip altitude, nautical miles 

Figure 1.- Angle -of -at tack oscillations and static trajectory of a skip rocket having an L/D 

of 6 and an over-all range of 4000 nautical miles. 
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Range, feet x 10" 4 
(a) Static trajectory. 

Figure 2,- Effect of gravity on trajectory and. oscillatory history of a skip rocket; L/D = 6 

over-all range = 4000 nautical miles. 
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(c) Envelope of angle- of -attack oscillations 
Figure 2,- Continued. 
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(d) Oscillation frequency 
Figure 2.- Concluded. 
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(a) Static trajectory. 

Figure 3 .- Effect of gravity on trajectory and oscillatory history of a 
skip rocket; l/D = 3; over-all range = 4000 nautical miles. 
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(d) Oscillation frequency. 
Figure 3 • ~ Concluded. 
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